Boltzmann equation for dissipative gases in homogeneous states with nonlinear friction 
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Combining analytical and numerical methods, we study within the framework of the homogeneous 
non-linear Boltzmann equation, a broad class of models relevant for the dynamics of dissipative flu- 
ids, including granular gases. We use the new method presented in a previous paper [J. Stat. Phys. 
124, 549 (2006)] and extend our results to a different heating mechanism, namely a deterministic 
non-linear friction force. We derive analytically the high energy tail of the velocity distribution 
and compare the theoretical predictions with high precision numerical simulations. Stretched expo- 
nential forms are obtained when the non-equilibrium steady state is stable. We derive sub-leading 
corrections and emphasize their relevance. In marginal stability cases, power-law behaviors arise, 
with exponents obtained as the roots of transcendental equations. We also consider some simple 
BGK (Bhatnagar, Gross, Krook) models, driven by similar heating devices, to test the robustness 
of our predictions. 
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I. INTRODUCTION 

Granular materials represent one of the most well- 
known paradigms for open dissipative systems; their 
ubiquitous character in natural phenomena or industrial 
processes, the possibility to study them at both very ap- 
plied and very fundamental levels have prompted many 
efforts to understand their properties ]E S B B Hf, 
but they remain challenging from many points of view. 
Thermodynamic-like descriptions remain in particular 
elusive for these intrinsically far from equilibrium sys- 
tems, because energy is continuously lost through in- 
ternal dissipation, and has to be compensated by non- 
thermal sources. 

Dilute granular gases present a particularly interesting 
framework which can be studied by model experiments, 
numerical simulations, hydrodynamics, kinetic theory or 
more phcnomenological approaches. The investigation of 
the velocity distribution of the particles reveals a rich 
phenomenology, with strong deviations from the equi- 
librium Maxwell-Boltzmann behavior, and its descrip- 
tion is still an object of debates. From an experimen- 
tal point of view, the situation may appear confusing 
at first sight. Although the measured velocity distri- 
bution F(v) generically deviates from the Maxwellian, 
its functional form depends on material property and 
on the forcing mechanism used to compensate for col- 
lisional loss of energy @, 0, H, [§]• A similar picture 
emerges from numerical simulations and analytical stud- 




where in 

addition to the more common stretched exponential be- 
havior, power-law distributions have also been reported 

m m m, m, m . 

Since the universality of the Maxwell-Boltzmann dis- 
tribution for equilibrium gases seems to have no analog in 



steady states of dissipative gases, it is necessary to iden- 
tify the generic trends for the velocity distribution and 
study their properties in a unified and simple framework. 
Our purpose is to develop further an existing quantitative 
cartoon to unveil the different effects at work, that lead to 
the wide range of behaviors alluded to above, when cer- 
tain key parameters are changed. To this aim, we focus 
on the framework of the generalized inelastic Boltzmann 
equation, which describes dilute homogeneous dissipative 
gases. The study of homogeneous systems is not only a 
useful starting point ,it is also relevant to experiments 
with bulk driving 0, Q . Spatial heterogeneity (and thus 
gravity, hydrodynamic instabilities, shock formation and 
clustering Q) will henceforth be discarded. 

The most characteristic features of the velocity distri- 
bution in dissipative systems - whether observed in time 
dependent scaling states or in non-equilibrium steady 
states - is overpopulated high energy tails. Generically 
these tails are stretched exponentials [HI, 1£ , HH . The 



more spectacular power law tails [H, l20U2ll 22L l3l| are 
exceptional. They were only found in systems of Maxwell 
molecules, and in the unusual device of Ref.[31[. In gen- 
eral these power law tails do not evolve naturally, but 
require careful fine-tuning of the physical parameters, 
both in the interactions, and in the driving mechanisms 

mm 

In the present paper, specific features and properties 
of F(v) will be derived and analyzed, when energy is in- 
jected by a negative friction thermostat, parametrized by 
an exponent 9. An important new result of general im- 
portance (321 . 133] ] is also the direct and simple relation 
between the parameters controlling the stability of the 
energy balance equation, i.e. the balance between colli- 
sional dissipation and energy injection on the one hand, 
and on the other hand the occurrence of new power law 
tails - l/v a ^ in the velocity distribution of the non 
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equilibrium steady state, which ap pea r at the margin of 
stability in this equation. In Ref. [331 ] a short summary 
of our results was presented. In this longer paper we 
show how all these results, and in particular these new 
power law tails, have been obtained with the help of a 
new method developed in Ref. [12, , which is applica- 
ble to different types of forcing mechanisms, and to a 
large class of interaction models, including inelastic hard 
spheres and Maxwell molecules. It seems virtually im- 
possible to obtain the present results by using the old 
methods developed for Maxwell molecules [H, [2(| HI] • 

The current paper hence constitutes a sequel to our 
first study 32], in which we have developed a method for 
analyzing the deviations of F(v) from Gaussian behavior. 
While we have focused in (32j on granular gases for which 
energy is injected by random forces, we consider here a 
different driving mechanism. We show that the large 
velocity tail of F(v) is characterized by a stretched ex- 
ponential exp[— v b ] with an exponent 6, that governs the 
stability of the non-equilibrium steady state. When this 
state is a stable fixed point of the dynamics, the expo- 
nent satisfies b > 0. For the exceptional cases of marginal 
stability, where b vanishes, F is of power law type, and 
the corresponding exponents will be calculated. 

The paper is organized as follows: we recall in section 
HTl the Boltzmann equation for inelastic soft spheres, to- 
gether with the criterion of stability of non-equilibrium 
steady states (NESS). We also briefly recall in III Bl the 
method introduced in our previous paper [32j . Sections 
IIIII and IIVI give the results for the large velocity tail 
of the velocity distribution in the case of energy injec- 
tion by non-linear negative friction, for both stable and 
marginally stable NESS. We turn in section [V] to a sim- 
ple linear model which mimics the basics of the inelastic 
Boltzmann equation. While this BGK model is amenable 
to analytical treatment, we show that it does not repro- 
duce the rich behavior of the non-linear Boltzmann equa- 
tion, in particular it fails for hard interactions. Section 
IVII finally contains our conclusions. 

II. INELASTIC BOLTZMANN EQUATION 

We consider the Boltzmann equation for a gas of in- 
elastic soft spheres [32| , which represents one of the sim- 
plest models for rapid granular flows. The dynamics is 
described as a succession of uncorrelated inelastic binary 
collisions, modelled by soft spheres with a collision fre- 
quency and a coefficient of normal restitution a, where 
< a < 1 0- The collision law (vi,v 2 ) — > (v^Vj) 
reads: 

vi = vi -p(g-n)n, v 2 = v 2 + p(g-n)n (1) 

where g = vi — v 2 , p = 1 — q = ^(1 + a) and n 
is a unit vector parallel to the impact direction con- 
necting particles 1 and 2. Inelastic collisions conserve 
mass and momentum, and dissipate kinetic energy at a 
rate cx ^(1 — a 2 ) = 2pq (the elastic case corresponds 



to a = 1). We consider a general collision frequency, 
<7?(<7,$) ~ 9 v \g- n | CT - Here ?(<7,i?) is the differential scat- 
tering cross section with ■& = cos _1 (g • n), v describes 
its energy dependence, and a its angular dependence. 
The exponent a ^ 1 describes a distribution of impact 
parameters biased towards grazing [a < 1) or head-on 
(<r > 1) collisions. For mathematical convenience mod- 
els with a = v have also been considered [2l|, |22j, |32j]. 
The symbol a = a/a denotes a unit vector, parallel to a. 
For elastic particles interacting via a soft sphere potential 
U(r) oc r~ a , one has v = 1 — 2(d— l)/a [34], where d is the 
space dimension. The exponents v = a = 1 correspond 
to standard hard-sphere behavior (a — > oo), while v = 
corresponds to Maxwell molecules (a = 2(d — 1)). Here 
v and a will be free exponents, that parametrize the ma- 
terial properties together with the inelasticity parameter 
a. 

We now give a simple representation of the nonlin- 
ear Boltzmann collision operator, which is convenient 
to study the spectral properties of the linearized colli- 
sion operator. The time dependent distribution -F(v, t) 
in spatially homogeneous systems obeys the nonlinear 
Boltzmann equation [331 ]: 

d t F{v,t)+TF = I{v\F) = J n J dvidv 2 ^|g-n| CT x 

F( Vl , t)F(v 2 , t) [5(v - v[) - <J(v - vx )] (2) 

where the collision operator I(v\F) has the usual gain- 
loss structure. For anisotropic F(v,t) the angular in- 
tegral, /„(■■■) = J*' 'dn(---)/ ' dn denotes an av- 
erage over the pre-collision hemisphere, g • n < 0. For 
isotropic distributions, as considered here, the integrals 
over pre- and post-collision hemisphere are the same, and 
/ (• • • ) = / dn(- • • )/fld can be extended over the com- 
plete solid angle, i.e. where = 2ir d / 2 /T{d/2) is the 
surface area of a d-dimensional unit sphere, and T(x) is 
the Gamma function. 

The forcing term TF represents the energy supply, 
working against inelastic dissipation. It may lead to a 
NESS (Non- Equilibrium Steady State). Absence of forc- 
ing {T = 0) describes free cooling, where the energy is 
decreasing in time. A heating device, considered fre- 
quently, consists in a random force acting on the par- 
ticles in between collisions: the corresponding stochastic 
White Noise (WN), widely used in analytical and numer- 
ical studies [SEalillMEiii, is described by adding 
a diffusion term —Dd^F to the Boltzmann equation. Our 
previous paper [32| has focused on this case. An interest- 
ing alternative is given by deterministic nonlinear Nega- 
tive Friction (NF). In general, the forcing term reads 

TF = (^ NF + J^wn)^ = 7<9v • (vv 6 -^) - D d?F (3) 

with (7 = 0, D ^ 0) for WN and (7 > 0, D = 0) for NF, 
and 9 > 0. The value 9 = 1 corresponding to the Gaus- 
sian thermostat allows us to study the long time scaling 
regime of an unforced system (so-called free cooling) (see 
section III Al and [iH]). The special case 9 = models 
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gliding friction. In general 9 is a continuous exponent se- 
lectively controlling the energy injection mechanism. In 
summary, the 'phase space' to explore is thus given by 
the parameters [y, a, 9, a) for NF driving, which is the 
main focus of this paper. 

The explicit form of I(v\F) above is conve- 
nient for calculating the rate of change of averages, 
/ dv^(v)F(v,t) = (tp\F) t , i.e. 



dt (1>(y)\F)t + (^F) t = \ J J d Vl dv 2 g»\g • nr x 

F(v 1 ,<)F(v 2 ,t)[V'(vi)+V(vi)-^(v 1 )-V(v 2 )] (4) 

Here ([2]) is just a special case of (HJ), where (ip(v)\F) t = 
F(w,t) for ip(v) = 5(v — w). The loss rate of energy 
follows by setting (tp(v)\F) t = {\mv 2 \F) t = |mu§(t), 
where Vo(t) is the r.m.s. velocity and v 2 the granu- 
lar temperature. We further use the normalizations, 
((l,vy)\F) t = (l,0,dv 2 (t)/2). 

A. Scaling form and stability of steady states 

To analyze the behavior of the distribution function we 
assume a rapid approach to a scaling form, 

F(v ) t) = (l/v (t)) d f(v/v (t)) (5) 

(for related proofs, see [24 ]) with normalizations 

<(l,c 2 )|/) = (l,d/2). (6) 

The inelastic Boltzmann equation @ can then be de- 
coupled into a time-independent equation for the scaling 
form /(c) and a time-dependent equation for the r.m.s. 
velocity vq, or granular temperature Vq, which reads 



d ( d 



dt \2 



(7) 



The collisional average (v 2 \I)t follows from (j4|) and ([T]) 
by carrying out the angular average, and by inserting ([3]) 
to change to scaling variables. Here J |g • n| CT+2 = /3 CT+2 
is given by 

& = / |g • n\° = r(2±i)r(|)/r(^)r(i) (8) 

J n 

where (3 a is well defined for a > — 1. The forcing terms, 
{v 2 \TF) t are calculated along the same lines using © 
and partial integrations. The results are, 

(v 2 \l)t = -±\2{(g» +2 )K +2 (t) 



(v 2 \F WN F) t = -2dD 
(v 2 \^ F F) t = -2 7 (c e+1 )« e+1 (i) 



(9) 



where g = \c\ — Ca|. The coefficient A 2 = 2pq(3 a+ 2 can 
be identified as the eigenvalue of the linearized Boltz- 
mann collision operator (see the general expression for 



A s below, Eq. ( (|18[) ) next subsection). Here, the no- 
tation ((fc(c!,c 2 ))) stands for the average of a function 
fc(ci,c 2 ) with weight /(ci)/(c 2 ). 

Let us first consider forced systems. The two terms in 
([7]) can then balance each other and lead to a NESS. For 
the WN-driven case (7 = 0), Eqs.flU)-© give 

lA 2 ((^+ 2 )X+ 2 (i) 



dv Q J dt 



4D 



AD 



1 



vg(t) 
■uo(oo) 



2// 



(10) 



where 6wn = 1 + an d vo(°°) is defined as the station- 
ary solution of PU|) . Similarly we obtain for the NF-case, 



dvo/dt = (^)( c «+i)^(t)-^ (( ^ 



1 



vo(t) 

Do(oo) 



)K 

r 



(t) 



(11) 



where &nf = v + 1 — 9. The dynamics always admits 
a fixed point solution of the equations above. The fixed 
point solution «o(oo) is stable/attracting for b > 0, unsta- 
ble/repelling for b < 0, and marginally stable for b = 0. 
As long as b > 0, the system naturally evolves towards 
the stable NESS. Note that D and 7 in ((lQj) and dTTJ) 
are irrelevant phenomenological constants, that can be 
absorbed in energy and time scales. 

In the stable NESS (vq — 0, b > 0), the corresponding 
integral equations for /(c) follow from ((2|) with dtF = 
and ((3]), and yield respectively for WN and NF, 



I(c\f) 
Hc\f) 



' (t>o(oo)) 2i> 

-d c -(cc s f) 



d 2 c .f = -Mg" +2 ))d 2 f (12) 



T..„. -x ,i' 7 " c l " ' ' ~ A 4 < (i 9 «+i) >> ^c-(cc e /)(13) 

Here -D and 7 have been eliminated with the help of the 
steady state solutions of (flC)]) and (fTTj) . 

It is also possible to consider the freely evolving state 
(FC), that does not lead to a NESS since 7 = and 
D = 0, but to a scaling solution /(c). Here the r.m.s. 
velocity decays according to 

*o = -t((9' /+2 )K +1 (t). (14) 

Moreover comparison of the FC case with the Gaussian 
thermostat (NF: 9 — 1) shows that the correspondin g in - 
tegral equations for /(c) are identical, as observed in [16j . 
This occurs because the term arising from dtF in ([2]) for 
FC is non-zero (since we do not have a NESS) and cor- 
responds exactly with the forcing term in Eq. Q13p when 
9=1 (since (c 2 ) = d/2). The Gaussian thermostat in 
fact describes a system driven by a linear (negative) fric- 
tion force, a = 7 c. This corresponds to a linear rescaling 
of the velocities, which leaves /(c) invariant. While the 
NF with theta = 1 and FC states are equivalent at the 
level of the scaled velocity distribution function, they dif- 
fer in the evolution equation for the temperature. It is 
therefore not possible to extend the stability criterion of 
NF to FC systems. 

We also note that (jT4j) gives a generalization of Haff 's 
law, that describes the decay of the energy or granular 
temperature for inelastic soft spheres. For instance, for 
v > one finds v 2 {t) - t- 2 / y (for more details, see [32j). 



B. Linearized Boltzmann operator 

We now recall briefly the method introduced in 32], 
which allows us to obtain the behavior of the velocity 
distributions at large velocities. Suppose that for large 
c = v/vq the velocity distribution can be separated into 
two parts, /(c) = fo(c) + h(c), where h{c) is the singu- 
lar tail part that we want to determine, and /o(c) the 
presumably regular bulk part. The tail part h(c) may 
be exponentially bound ~ exp[— (3c b ] with < b < 2, or 
of power law type. In the bulk part /o(c) the variable 
c is effectively restricted to bulk values in the thermal 
range v < vq or c < 1. As far as large velocities are con- 
cerned, the thermal range of / may be viewed to zeroth 
approximation as a Dirac delta function 5(c), carrying all 
the mass of the distribution. In this way, we obtain an 
asymptotic expansion of /(c) by considering the ansatz 

f(c) = S(c) + h(c), (15) 

and linearizing the collision term I(c\f) around the delta 
function (using the relation I(c\d) — 0). This defines the 
linearized Boltzmann collision operator, 

I(c\S + h) = -Ah(c) + 0(h 2 ). (16) 

Note that we restrict the analysis to isotropic functions 
h(c). The eigenf unctions of A decay like powers c~ s . 
Consequently they are very suitable for describing power 
law tails, /(c) ~ c~ s . We note, however, that the mo- 
ments J dcc s f(c) are largely determined by the regular 
bulk part, which can therefore not be approximated by 
<5(c) when calculating moments. 

The most important spectral properties are the eigen- 
values and right and left eigenfunctions, 

Kc- s - d - v = X s c 8 "' ; AV = \ s c s+u (17) 

with eigenvalues for s > [12], 

A s = p a {1 -2-Fi (-§ , 2±L-^\l- q 2 )}- p s B s+a . 

(18) 

Here 2F\(a,b]c\z) is a hyper-geometric function and B a 
is given by l|8]). The value s = with Ao = is an 
isolated point of the spectrum with corresponding sta- 
tionary eigenfunctions, invariant under collisions, 

A<5(c) =0 and A f • 1 = (A = 0). (19) 

Right and left eigenfunctions are different because A is 
not self-adjoint. There is in fact among the eigenfunc- 
tions in (|17p another, less trivial, stationary right eigen- 
function, c~ s ~ d ~ u ; where s* is the root of transcendental 
equation A s = (see [Hj]). 

We also note that the eigenvalue X s is independent 
of the energy exponent v: it is the same for inelastic 
Maxwell molecules, hard spheres, very hard particles, 
and very weakly interacting particles. The reason is pre- 
sumably that the scattering laws are the same in all mod- 
els, and equal to those of inelastic hard spheres. More- 
over, it depends strongly on the inelasticity through a, 
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FIG. 1: Concave eigenvalue spectrum A s (cr) for s > of the 
collision operator in the inelastic soft sphere models {a, v}, 
defined in the text, and shown for various values of the co- 
efficient of restitution a. The ordinate shows A s (cr)//3 CT for 
d = 2, a = v = 1, which approaches f for s — > oo, and -1 
for s — * 0. The point Ao(cr) = is an isolated point of the 
spectrum. 

and weakly on the angular exponent a. Fig. [T] shows 
that A s is a concave function of s. 

The previous paragraphs refer to the action of the 
I{c\f) on functions /(c) of power law type. We shall also 
need the large-c form, 1^, of the collision operator act- 
ing on exponentially bound functions, /(c) ~ exp[— /?c fc ], 
with positive constants /?, b. The nonlinear operator re- 
duces to a linear one for c\ 3> Ca, because g ~ ci and 
c[ = c —p(ci ■ n)n, and the C2 -integration can be carried 
out in ([5]). The resulting expression is, 

4o(c|/) = -A,c"[l - /C CT c- b ( CT+1 >/ 2 ]/(c), (20) 

and we obtain for the coefficient, 

K a - (r(^)/r(^)) ((2//36(l - g 2 ))^ 1 )/ 2 ) (21) 

as derived in Appendix B of Ref. [13] ■ Note also that this 
coefficient vanishes in one dimension. 

III. ASYMPTOTICS FOR STABLE NESS 

Exact closed forms of the scaling solution of (|f 3[) are 
not known in general, but the high energy tail may be 
computed accurately with the method developed in Ref. 
f32l ]. Restricting the analysis to standard arguments 
where the asymptotic form of I(c\f) is Joo(c|.f) ~ — B^c" 
(see e.g. [HI), an interesting feature emerges [32|, [33|: in 
the region of stability (b > 0), the asymptotic solution 
of (fT3|) has a stretched exponential form, / ~ exp[— c b ], 
with b = v + f — 9, while in cases of marginal stability 
(6 — * + ), / ~ cr a is of power law type with an a priori 
unknown exponent a. As expected, b decreases when v 
decreases, since a tail particle with velocity c 3> f suffers 
collisions at a rate c v . The slower the rate, the slower the 
particle redistributes its energy over the thermal range 
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c < 1, which results in an increasingly overpopulated 
high energy tail. When v is further decreased such that 
b changes sign, the tail is no longer able to equilibrate 
with the thermal "bulk" , and the system cannot sustain 
a steady state. A similar intuitive picture may be devel- 
oped with respect to 9 in the NF cases (33J. 

The leading behavior of /(c) is a generalization of 
the one obtained for inelastic hard spheres (v = 1) and 
Maxwell molecules (y = 0). While we had analyzed in 
details the case of WN driving in [33] , we focus here on 
inelastic gases driven by NF, which include the Gaus- 
sian thermostat, or equivalently, the Freely Cooling gas 
(FC: 6 = 1). The method allows us to calculate the 
sub- leading correction for c > 1 to Ioc(c\f). This in 
turn yields important sub-leading multiplicative correc- 
tion factors to /(c) of exponential and power law type, 



In /(c) ~ -(3c b + f3'c b ' + X In c + 0(1), 



(22) 



where b > b' > 0. This expression is in fact an asymptotic 
expansion of In /(c). The limiting corrections as b' — > + 
are already contained in the exponent %. Moreover, as 
soon as b' becomes negative, the correction term becomes 



<C C(l), and should be neglected for consistency. 



In the spirit of asymptotic expansions we only look for the 
sub-dominant correction c b with b' > and set x = 0. 
Only if b' = do we look for terms of type x m c - The 
goal of this section is to calculate the exponents {b, b', x} 
explicitly, and to express the coefficients /3'} in terms 
of the moments ((g v+2 )) and (c 6+1 ). These moments can 
be independently measured in the DSMC (Direct Simu- 
lation Monte Carlo) method (see Ref. [32|]). The sub- 
leading approximation supposedly extends the agreement 
of theoretical predictions with measured DSMC data to 
smaller c- values. 

We start with the NF integral equation, obtained from 
(fTB")) by replacing the collision operator / by the full 
asymptotic form 1^ in (|20p . The last form is the appro- 
priate one for exponentially bound functions /(c). This 
yields 

-B^ii-ic^c-^ + ^ b }f = ^f+id+e-iy-'f, 

(23) 

where the c-independent factors have been combined 
into, 



2(d + a)(c e 



\ 2 ((g»+ 2 )) (l + a)pq((9»+ 2 )y 



(24) 



and A2 = 2pqf3 a+2 and (J8]) have been used. Here the con- 
stant B a depends on all three model parameters (y, a, 9), 
and contains averages with the unknown weight /(c). 

The parameters in /(c) can be obtained from the full 
integral equation ([23]) by substituting the ansatz (|22)l . 
applying the derivative, equating leading and sub-leading 
powers of c, and recalling the relation b > b' > 0. To 
leading order we have 6/3c f,+9_1 — c v B a , yielding 



1 



(3b = B a 



(25) 



The exponent b is the same as the one found in the sta- 
bility analysis in (fTTj) . These results are largely gener- 
alizations of special cases, existing in the literature for 
9 = {0, 1}; v = {0, 1}; a = {1, v}, derived in 0. 

The remaining terms with sub- leading powers of c have 
respectively the exponents E\ = b' + 6 — 1,E$ = 9 — 
1, E3 = v— i(cr+l)fo. First consider the case (7 = 1, where 
£3 = v - b = 9 - 1 = E 2 . As Ei > E 2 the coefficient 
(3' = 0, equating the coefficients of the remaining terms 
then yields the second line of the equation (f2T)|) below. 
If er > 1, then E\ > E 2 > E3, and the coefficient of 
each power has to vanish, yielding the first line below. If 
er < 1 we obtain the sub-leading term by matching the 
exponents E\ = E3, yielding the third line below. Lower 
order terms with exponents E 2 have to be neglected for 
consistency, hence x = 0. So, the sub-leading results for 
NF driving are, 



a > 1 

(7=1 
(7 < 1 



13' = 0, x=l~d-0 
/3' = 0, x = l-d-9+(3b)C 1 (d) = -9 



(d-l)q 2 



0, b' = ±b(l-a), f3'b' 



l-q z 
f3bJC a {d). (26) 



In one dimension the results simplify substantially. The 
collision kernel in the Boltzmann equation @ lacks the 
angular integration J , fi a = 1 in ©, and K-i(d) = in 
([2T|l for all cr, implying f)' = 0, and b' is irrelevant. Then 
the large-c behavior of the distribution function is, 



f(c) 



'exp[-/3c 6 ] (d=l). 



(27) 



Other simplifications occur for special values of the pa- 
rameters 9 and v. For v = (Maxwell molecules) the 
coefficient in JM]) simplifies as ((g l/+2 )) = ((g 2 )) = d. 

Simplification also occur for a case of special interest, 
the free cooling system or equivalently the Gaussian ther- 
mostat (9 — 1), where (c e+1 ) = (c 2 ) — ^d on account of 
([6]). Here the exponents and coefficients are to leading 
order (see ([Ml)). 



b = v, f3b = B a 



and in sub- leading order, 



(d + a)d 



(l + o-)pq((g»+ 2 )) 



(28) 



cr > 1 : /?' = 0, X=-d 

dq 2 - 1 



(7=1: p' = 0, X = 



l-q* 



(29) 



<7<1: x = 0, b' = ifo(l-(7), 0'b' = PbtCid). 

For free cooling [9 = 1) at d = 1 we have /Ccr(l) = 0, 
hence (3' = 0, yielding, 



/(c) ~ (l/c)exp[-/3c fc 



(30) 



Further simplification occurs in freely cooling Maxwell 
models {v = 0, 9 = 1), where B a = (d + cr)/[(l + a)pq ]. 
This is a marginally stable case (b = v+ l — 9 = 0), and 
will be discussed in the next section. 
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Finally we compare the analytic predictions with the 
DSMC results. The DSMC method offers a particu- 
larly efficient algorithm to solve the nonlinear Boltzmann 
equation [36| . Figure [5] shows for the one-dimensional 
case the simulation results (solid line) for free cooling 
(9 = 1) in the soft sphere model [y = i) with com- 
pletely inelastic collisions (f> = q = i), compared with 
the analytic results in zeroth approximation (dashed- 
dotted line), i.e. / ~ exp[— f3^/c\ in ([28]) . and in first 
approximation (dashed line), cf ~ exp[— (3^fc\ in (|30|) . 
where — 8/((g 5 / 2 )) (according to ([22])) is obtained by 
an independent DSMC measurement of the two-particle 
moment. The zeroth approximation has an effective slope 
different from the slope of the first approximation. The 
latter essentially coincides with the DSMC measurements 
for all c > 1.7. Also note that the theoretical curves can 
be shifted in the vertical direction to give the best pos- 
sible fit with the DSMC data, because the overall con- 
stant factor exp[C(l)] in /(c) cannot be determined in 
oui 
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FIG. 3: Free cooling with (0 = l,d = 2,6 = v = a = 0.5) at 
various values of a. 



enough. However, the solid curve (transformed DSMC 
data) /(c) exp[— (3'c b ] shows a striking agreement with 
the theory exp(— (3c b ), and demonstrates that the sub- 
leading corrections extend the validity of the asymptotic 
theory to much smaller c— values, thus enabling us to 
test the validity of theory, and establish the importance 
of the sub-leading corrections. Striking is the fact that 
such plots of /(c) vs c h produces linear high energy tails 
(in spite of the importance of the sub- leading correction) , 
which would then be well fitted with an effective value of 
[3: /(c) ~ exp(— (3 c sc b ) (this is also the case in Fig. [2]). 
As shown here, such an effective value can be markedly 
different from the true /?, which indicates that any fitting 
procedure, aiming at computing 0, is doomed to fail. 



FIG. 2: Free cooling with (0 = 1, d = 1, b = v = 1/2, a = 0), 
where /(c) vs. c b and c/(c) vs. c 6 , are compared with 
exp( — Pc b ), to show the exp( — f3c b )/c behavior of /. The 
solid line represents the Monte Carlo (DSMC) data. The 
inset shows the overpopulation of the high energy tail when 
compared to a Gaussian (on such a plot, a Gaussian would 
produce a concave instead of convex graph). 

The DSMC data in Fig. [3] at large velocities show 
the stretched Gaussian behavior exp[— 0y/c] for two di- 
mensional free cooling in the soft sphere model with 
(b = a = v = i). They indicate that the coefficient (3 in- 
creases with a. This figure illustrates the overpopulation 
of the tail with respect to a Gaussian (indistinguishable 
from the a = 0.9 curve here, shown with stars). 

Striking examples of the importance of sub-leading 
corrections are shown in Fig. 21 for a two-dimensional 
model with 9 = 1, is = 2, a = and a = —0.5. In 
these cases b = 2, b' = 1 (for a = 0) and b' = 3/2 
(for a = -0.5). Comparison of the "raw" DSMC 
data (dashed curve) with the dominant asymptotic pre- 
diction exp(— (3c b ) (dotted curve) shows no agreement. 
The reason is that the simulated c— values are not large 
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FIG. 4: Comparison of the velocity distribution obtained from 
Monte Carlo (DSMC) simulations with the asymptotic predic- 
tions. For the Gaussian thermostat (6 — 1, d = 2, a = 0, v = 
2) at a — (top) and a — —0.5 (bottom) the exponents are 
(b,b', X ) = (2,1,0) and (2,1.5,0) respectively. DSMC data 
are plotted as /(c) (dotted line) and exp[— (3'c b ]f(c) (solid 
line) vs x = f3c b , and compared with the theoretical predic- 
tion e~ x (dashed line). Here (j3, /?') ~ (1.087, 1.359) for a = 0, 
and ~ (1.585, 1.616) for a — —0.5 have been measured in the 
DSMC simulations from their definition given in the text. 
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Figure [5] shows DSMC data for NF forcing with var- 
ious values of v and 9. The simulations confirm the 
large-c predictions, i.e. In /(c) ~ —(3c b , where b = 
v + 1 — 9. Moreover, the dashed lines show the agree- 
ment with the prediction exp(— (3c b ) where the coefficient 
(3 = 2(\c\ e+1 ) /[bpq{(g u+2 ))}. For those parameters, the 
sub-leading corrections are negligible. 




— v=0 
v=0.5 

— v=l 

— - v =2 



10 



15 




FIG. 5: Negative friction with (6 = v + 1 - 0, d = 1, a = 0). 
Plots show /(c) vs. c b (A) for various values of v at 6 = 
0.5, and (B) for various values of at v — 0.5. The dashed 
lines correspond to the predictions exp(— /3c 6 ) vs. c b with (3 
calculated in each case by DSMC. The inset corresponds to a 
two-dimensional case {6 = 0, v = 2, a = 0), showing /(c) vs. 



Regarding the soft sphere systems in stable NESS 
(b > 0), either freely cooling or driven by Gaussian ther- 
mostats, we may conclude that the agreement between 
analytic and DSMC results for high energy tails is very 
good. 



MARGINAL STABILITY, POWER LAW 
TAILS 



IV. 



We now analyze the integral equation (|T3|) for the 
threshold models (b — 0; this fixes the exponent v at 
the threshold). Marginal stability is a limiting property 
of a stable NESS as b -> 0+, which occurs in states, 



driven either by white noise (see Ref.[32]) or by negative 
friction. 

As we have seen in the previous section, the high en- 
ergy tails for stable states (b > 0) have the generic form 
/(c) ~ exp[— /3c b ] with j3 = B a /b, and sub-leading correc- 
tion factors of similar structure. To illustrate how power 
law tails arise, we take the limit of /(c) as b — > + using 
the relation (c b — l)/b ~ Inc. The result is, 



f (c) ~ lim a 

1 ■ ' 6^0+ 



exp[-£U76] 



(31) 



Of course (\ — rf) is not the full exponent of the tail, 
because the exponential form above represents only the 
leading asymptotic behavior for b > 0. For instance, any 
correction factor exp[— j3'c h ], where b' = H{b) — > as 
6^0, gives additional contributions to (pTTj) . 



A. Gaussian thermostat (NF: 9 = 1; b = u = 0) 

Here the Maxwell molecules are the marginally stable 
model. To determine the full exponent of the power law 
tail we linearize the nonlinear integral equation (|13p at 
the stability threshold around the "thermal bulk part" 
of /(c), using (Tr5|) and (TTo| . We start with the simplest 
case of inelastic soft spheres, driven by a linear friction 
force (9 = 1). 

Substitution of /(c) = 6(c) + h(c) in the collision kernel 
of (fT3"|) yields to linear order in h(c), I(c\S + h) = —Ah(c). 
The r.h.s. of (H]) also simplifies, as {{g v+2 )) = {{g 2 )) = 
d, and the resulting integral equation is, 



Ah 



(c/)- 



(32) 



Inspection of this equation and (fT7|) shows that the op- 
erators on left and right hand side, when acting on the 
right eigenfunction l/c s+d (recall that v = 0) generate 
new powers of c. Solving the integral equation implies 
that one determines the value s* that makes both expo- 
nents equal, leading to the transcendental equation, 



(33) 



Consequently, the solution of (f3"2"| . which presents the 
asymptotic large-c solution of (|13p , is the power law tail, 



/(c) ~ h(c) ~ I/c 



s*+d 



(c»I). 



(34) 



If the transcendental equation has more solutions, then 
the largest root s* is the relevant one, because the energy 
(c 2 ), and possible moments (c°) and ((g a )}, appearing in 
the transcendental equations (see next subsection) must 
be finite, imposing s* > max{2,a}. So, the obvious so- 
lution of (|33| . s* — 2, has to be rejected. However the 
equation has a second solution with s* > 2, because A s 
is a concave function of s. This can be seen directly 
from a graphical solution by adding in Fig. [T] the line 
y(s) = 7jsA 2 . The numerical values of s*(a), obtained 



from the numerical solution of (|33j) , are shown in the in- 
set of Fig. for the two-dimensional system. The main 




FIG. 6: Power law tails in free cooling, obtained for the 
threshold model (6 = 1, d = 2, a = 1, b = v = 0). The inset 
compares predicted and measured exponents. As a increases 
the exponent increases, and the curve tends to a Gaussian. 



plot shows the comparison of the DSMC measurements 
of /(c) for this system compared to the theoretical pre- 
dictions. 

It is also instructive to consider the one-dimensional 
version of ([55]) . which can be solved analytically. Then 
the eigenvalue (fT8|) simplifies to A s = 1 — q s — p s , and (|33|) 
becomes, 1 — q s — p s — spq, with solutions s* = {2, 3}, and 
s* = 3 is the relevant one, and /(c) ~ l/c s +d ~ 1/c 4 in 



_2T2 



agreement with the exact solution /(c) = (2/7r)/[l + c 
found in [l3l |. 

For <i > 1 the transcendental equation can not be 
solved analytically, except in a few limiting cases, that 
we consider first. In the elastic limit (a — ► 1 or q — » 0) 
one only finds a meaningful solution of (|33() by letting 
simultaneously s — > oo while keeping sq = £ = fixed. As 
pr+z/p" = (cr + l)/(cr + d), and A S //3 CT -> 1 as s -> oo 
(see Fig. Q]or Eq.(3.12) in Ref.JIl]), the transcendental 
equation (|4Tj) approaches 1 ~ £(1 + cr)/(d + a), yielding 
the solution, 

8%=t/q~[(d + a)/(l + o)]/q (a -> 1) (35) 

for general cr. In the elastic limit as a = 1 — 2q — > 1, the 
root s* moves to oo and the algebraic tail disappears, as 
required by consistency with the Maxwell distribution in 
the elastic limit. Using the large s— expansion of (|33| it 
is straightforward to obtain additional sub-leading cor- 
rections. 

Another case where the integral equati on (l33l) can be 
solved analytically is at large dimensions [2l[. To do so 
it is convenient to divide (f3"3"| by j3 a . As d — * oo its right 
hand side approaches spq(l + a)/d. So, one finds only 
a meaningful solution by simultaneously letting s — > oo 
while keeping x = s/d = fixed. To calculate \ s / Pa from 
(|18[) in this coupled limit we use the relation, 

lim 2 F x (-f , *±i; 2±S|z) = E:=o(^(-^)' 
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FIG. 7: (a, d) /d vs a for the two variants of d~ dimensional 
Maxwell models, i.e. {top) with (a = 1, v = 0) and (bottom) 
with (cr = i/ = 0). The dashed lines correspond to the analytic 
results (|38[) for large d. 



where (a) n = T(a + ri)/T(a). This relation can be de- 
rived starting from the Gauss hyper-geometric series (3?| 
for 2-Fi(a, b; c\z) by taking the (d — * oo) limit term by 
term, and subsequently using the relation 2^1 (a, b; b\z) = 
(1 — z)~ a . Then (f33|) for the present threshold model sim- 
plifies to, 



1- (1 + 3(1 -q 2 )) 



-(l+a)/2 



xpq(l 



(37) 



For the cr-values, mostly considered in the literature, i.e 
1) " 



the model with [a = 1) 2 
convenient model (a = 0) 



and the mathematically 
22j, the above equation 



can be solved analytically. For the Maxwell model with 
a = 1 it is a quadratic equation, and for a — it is a 
cubic equation. The root x = is not a solution of (f33|) 
because X s in (fl8| holds only for s > 0. The resulting 
s* in (|34jl becomes in the coupled limit d — > 00, s — » 00 
with s/d = x = fixed, 



= xtd 




VT+zx' 

(36) 



(a = l) 

(a = 0) ' 
(38) 

The exponent s* + rf, obtained here, disagrees with the 
result of Ref. [2l| in the sign in front of the square root. 
We note that the a— dependence of and s\ in the last 
equation is somewhat different at large d. The exponents 
in ([33]) and (f3"5|) agree in the respective limits d — ► 00 and 
a — * 1. 

Equation (|33[) can easily be solved numerically. For 
the Maxwell model with a = 1 the resulting exponents 
s* and Sq as a function of a for various d are plotted as 
s\/d and Sq/c? in Fig. [7] As shown in the inset of Fig. [51 
the agreement with DSMC simulations is very good. 

Most results of this subsection, applying to Maxwell 
models (y — 0), have been derived already in the lit- 
erature using an entirely different mention, namely by 
Fourier transformation with respect to the velocity vari- 
ables. The Fourier transform method can only be applied 
to Maxwell models (y = 0) where the microscopic colli- 
sion frequency is independent of the relative velocity g, 
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leading to a collision kernel I(c\f) that is a convolution 
product in velocity space. It simplifies to an ordinary 
product after Fourier transformation. The method can 
not be generalized to inelastic soft sphere models with 
v 7^ 0. Our method on the other hand can be applied for 
all values of v. 



B. Nonlinear negative friction 
(NF: 8>Q;b = v + l- 6 = 0) 

In this case, the threshold model is the soft sphere 
mode with b = 0ovv = 9 — \. The corresponding scaling 
equation for the high energy tail is obtained by setting 
v = 9 — 1 in (|13p. and reads, 

I(c\S + h) = -Ah = \\ 2 T{6)d ■ (cc e h), (39) 

where we have defined the ratio of the moments T(6) as, 

((q e+1 )) 

m s (40) 

This quantity should not be confused with the Euler 
Gamma function. We also note that T(9) is unknown 
a priori, as it depends on the full unknown scaling form 
/(c) with c <E (0, oo). Inspection of ([39]) shows again 
that the operators on left and right hand side of (|39[). 
when acting on the right eigenfunction \/c s+d+l ' with 
v = 6 — 1, will produce new powers of c, and one deter- 
mines the value s* , that makes both exponents equal, by 
solving the transcendental equation, 

X s = isA 2 r(0) = S pq(3 a+2 T{9). (41) 

We recall that A s is the same for all inelastic soft sphere 
models. We further note that T(9 = 1) = 1, as can be 
verified from (|40|) and the normalization (c 2 ) = d/2, and 
we recover the transcendental equation (|33p for linear 
friction. 

Denoting the relevant root of (|4"Tj) by s* the solution of 
((39)) is the right eigenfunction of A with eigenvalue A s « , 
i.e. 

f(c)~h(c)~c- s °- d - e+1 . (42) 

So at the stability threshold for driving by nonlinear fric- 
tion (y = 9 — 1), there exists again a power law tail in 
the scaling solution of the Boltzmann equation for soft 
sphere models, provided (|4"Tj) does indeed have a real pos- 
itive solution. 

Extracting the largest root from (|4ip is somewhat more 
complicated than in equation (|33[) . because of the un- 
known factor r(#). Even for d = 1 there are no simple 
exact solutions. To obtain T(9) we determine the mo- 
ments in (|4H|) and their ratio T(9) by direct DSMC mea- 
surements. The inset of Fig. [8] shows T(6) resulting from 
these measurements at a = in two dimensions. The 
plot shows that T(9) is an approximately linear func- 
tion increasing with 9. At this point, it is noteworthy 



that a Gaussian ansatz for the velocity distribution yields 
r(0) = 2^ s ~ 1 ^' 2 . This provides an excellent approxima- 
tion (not shown) (42)], which also coincides with the ex- 
act value at 9 = — 1. There the linear approximation 
is slightly off. However, in the physically relevant range, 
9 £ [0, 1], the linear approximation is slightly better. The 
following analytical results confirm this trend: T = i, 1 
for 9 = —1, 1 respectively. The resulting T(9) is used as 
a known input parameter in (|4ip . 

Once T(9) is known from DSMC measurements, one 
can construct a simple graphical method for solving (|41j) 
and classifying its possible solutions for different values of 
9 and a. Here we discuss only the v models with a = 1. 
This is done by plotting in Fig. [5] for a fixed value of a 
the curve, y%{s) — \ s /fli (1-h.s. of f4pj)). and the straight 
lines, y 2 (s) = [2pqT(9)/(d + l)]s (r.h.s. of gTJ), as fol- 
lows from P3/P1 = 2/(d + 1)), for different values of 9, 
and determine the largest intersection point. Here T c de- 
fines the slope of the line, y = [2pqT c / (d + l)]s, through 
the origin, that is tangent to curve A s . The largest in- 
tersection point of the eigenvalue curve with the line, 
labelled T(9 = 1) = 1, represents the graphical solution 
for the linear friction case, and the relevant root s* (a) 
has already been obtained in Fig. [H] for two dimensions 
[22I [40| . For the nonlinear case we obtain the following 
scenario. As 9 decreases from 1 to 0, the ratio T(9) de- 
creases from 1 to some value T(0) > 1/2, and the largest 
root s+ = s+(#, a) grows from s* to some value Sj_(0, a). 
As 9 grows larger than 1, T(9,a) increases from 1 to 
r c (a), and the largest root s*_ = s*_(9,a) decrease from 
s* to s*, as shown in Fig. [5] For T(9) > T c the root of 
the transcendental equation becomes complex. The cor- 
responding tail with an oscillatory pre-factor is no longer 
everywhere non-negative, and thus becomes unphysical. 

We finally present a comparison of our analytical pre- 
dictions with the result of DSMC simulations for sev- 
eral parameter values in Figs. [5] and [TO] Due to the 
marginally stable character of the NESS, simulations are 
quite difficult and time-consuming. Nevertheless, an ex- 
cellent agreement is obtained for all parameter values. 
Note that the values of the power-law exponents are large 
so that a direct fit to power-law forms would not be very 
precise, and could not exclude other fitting forms (since 
at most one decade in c is covered). 



V. INELASTIC BGK MODELS 

In this section we study a simple inelastic BGK (Bhat- 
nagar - Gross - Krook) model for homogeneous velocity 
relaxation [3^, H(|, which only takes the most essential 
features of the complex nonlinear collision operator into 
account. The goal is to understand how much of the 
rich behavior of the Boltzmann equation, described in 
the present paper and in [32j , is preserved in such a lin- 
ear model. The analytic results of the previous sections, 
and of Refs. [13, H3, HH], are restricted to asymptotic 
solutions, which can be applied directly to (14"3"1) . The 
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FIG. 8: NF at (d = 2, a = 1, u = 9 - 1, a = 0). Graphical 
solution of (I4ip for the marginally stable NF driven soft sphere 
model. The concave curve represents the eigenvalue yi(s) = 
A s (l)//?i (solid line) and the straight lines represent 1/2 (s) (see 
text) for different values of F(#), labelled from bottom to top 
by r n (n = 0, 1, 2). The inset shows F(9) versus 9 as obtained 
from DSMC measurements. The slope of the tangent line is 
labeled by F c . The largest intersection point corresponds for 
a given value of T{9) to the root a„(a), which determines the 
power law tail /(c) ~ l/c a+d+B - 1 . 




FIG. 9: NF at stability threshold (d = 1, 6 = v + 1 - 6 = 0). 
(top): (9,u) = (0.5,-0.5); {bottom): (9,v) = (0,-1). The 
lines are the predicted power law tails, /(c) ~ l/c a+d+0 ~ 1 , 
following from the construction discussed in Fig. E 

BGK kinetic equations allow to go further, since they 
reduce to simple linear first and second order inhomoge- 
neous ODE's, which can be solved exactly, at least for 
systems that are freely cooling, or equivalently driven 
by linear negative friction, as well as for systems driven 
by white noise. Although the present paper is mainly 
dealing with nonlinear negative friction, we restrict our- 
selves to the Gaussian thermostat (linear friction) and 
also discuss white noise driving for completeness. The 
exact solution of the BGK model with the full nonlinear 
friction is not known. 

In a crude scenario for relaxation without energy in- 
put, the velocity distribution F(v,t) relaxes towards 
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FIG. 10: NF at stability threshold (d = 2,b = v + 1 - 9 = 0). 
(top): a = and (0,v) = {(1, 0), (0.5, -0.5), (0.1, -0.9)}; 
(middle): (9,v) = (0.5,-0.5) and a = {0.0,0.3,0.5}; (bot- 
tom): (6,v) = (1.1,0.1) and a = {0.1,0.5,0.6}. The plots 
show the predicted power law tails, /(c) ~ l/c E , as dashed 
lines with exponents E = a + d + 9 — 1, where a is calculated 
from the transcendental equation (|41|l using the measured 
values of T(9) (see also Fig. 18} . The predicted exponents 
are: (top): E = {6.0, 6.8, 7.75} ; (middle): E = {6.8, 7.8, 9.4}; 
(bottom): E = {6.0,8.0,9.2}. These exponents show very 
good agreement with DSMC data. 

a Maxwellian with shrinking width avo(t), at a rate 
oc i>g(t). The width, proportional to a, mimics the role 
of the coefficient of restitution, which reduces the typical 
velocity in an inelastic collision by a factor a. With a 
constant supply of energy, the system can reach a NESS, 
and the global evolution can be modelled by the BGK- 
type kinetic equation, 

d t F(v,t) - Dd*F(v,t) = I(v\F) 

I(v\f) = -v v (t)[F(v,t) - F (v,t)\ 

F (v, t) = (av y d cf>{v/av (t)), (43) 

where (f>(c) = ir~ d / 2 exp[— c 2 ] is the Maxwellian. If F(v, t) 
is rapidly approaching the scaling form l[5]). the rescaled 
collision kernel, I(v\F) — v^ d I(c\f), takes the form, 

I(c\f) = -f(c) + a- d 4>(c/a). (44) 
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FIG. 11: BGK in one and two dimensions for FC: log-log plot 
showing the power-law tails. 



We note that the collision kernel does not show any v- 
dependence. This is similar to Maxwell models, where 
the collision frequency is independent of the microscopic 
velocities of the colliding particles. The time evolution 
of vo(t) in free cooling and in the WN (white noise) case 
obeys equations of motion, similar to ([?))- (fTTj) . Because 
for inelastic soft spheres the rhs is also proportional to 
1 — a 2 , the discussion about stability of the granular tem- 
perature T(t 
WN driving [3 



Vn(t) is the same as in free cooling and 



and the same applies to Haff 's ho- 
mogeneous cooling law, v§(£) ~ t~ 2 l v . 

In the free cooling case (D = 0) the energy equation be- 
comes, vq — — 2pqvQ +1 . Inserting then the scaling ansatz 
© in 62D yields, 



cf + {d+a)f=— d <t>(~) ] a -2/(1 -a 2 ) 
or \ol) 

Its exact solution is (see also [39l |40|). 

ft \ aa ° 
/(c) = 



(45) 



T d/2 



^d-\-a 



aa a T (^±2) 



2TT d / 2 



1 



(c » 1) 



(46) 



This solution, including its high energy tail (see Fig. ITT]) , 
is independent of the exponent v. A similar heavily over- 
populated power law tail, /(c) ~ l/c d+a with d > 1, is 
also found in the freely cooling Maxwell model. There the 
exponent a(a) takes in the elastic limit (a — ► 1) the very 
similar form a ~ 4d/ (1 — a 2 ) (compare (|45|l). We also 
note that the a-dependence of the power law exponent 
in the BGK model is essentially the same as for higher 
dimensional Maxwell models [32j , and similar to (|38|) for 
NF driving. However, in the general class of inelastic 
soft sphere models with collision frequency g v and v > 
(hard scatterers), the tail is not a heavily overpopulated 
power law tail, but a lightly overpopulated stretched ex- 
ponential, f(c) ~ exp[-/3c b ] with b = v > 0. The BGK 
models describe quite well the features of the soft scat- 
tering models, but are totally missing the more effective 
randomization caused by the high speed particles present, 
in models with positive v. 



Let us now turn to the case of white noise driving in 
the BGK model of Eqs. (|43[) . Again the energy balance 
equation is the same as (0 for soft spheres. So all BGK 
models with v > — 2 have a stable attracting fixed point 
wq(co), and the integral equation has a rescaled form, 
analogous to (Q]: 



/"(c) + — f( C )-2a/(c) = -^(^) 
c or \al 



2a 



(47) 



Here a prime on / denotes a derivative with respect to its 
argument c. Eq. (|37]l shows that /(c) is independent of 
the model parameter v and of the noise strength D. The 
equation can be solved exactly, and the two integration 
constants are fixed by the normalizations ([6]). For all 
values of d we make the transformation 



/(c) = a - d y(f3c); b = ±(3a; (3 = V2a~ = 2/Vl - a 2 , 

(48) 

where y is a function to be determined. The resulting 
equation for y can be solved: the one- dimensional BGK 
model has the exact solution 



i&exp[& 2 ] [e x eric{b+ §) + e^erfc^ - §)] 



(49) 

Using the properties of the complementary error function 
erfc(z) one can verify that the first term inside [• ■ • ] de- 
cays for x — * ±oo as exp[— a; 2 /46 2 ] and the second one as 
2exp[— x], yielding an exponential tail, 



/(c) ~ i/3exp[6 2 ]e 



(c»l). 



(50) 



Similarly we find in the two-dimensional case for the 
solution satisfying the normalizations ([B]), i-e. 



1 



K (x) / zdzexp[~z 2 /4b 2 ]I (z) 



1 f x 

+ -I (x) / zdzexp[~z 2 /4b 2 ]K a (z), (51) 
I" Jo 

where Io(x) and Kq(x) are Bessel functions with imagi- 
nary argument [13] • The exact solutions (|4*6")) , (|4*g|) and 
([STj) have been obtained by K. Shundyak (43|. At large x 
we have Ko(x) ~ e" x \JiT/2x and Io(x) ~ exp[— x 2 /4b 2 ], 
yielding the high energy tail, 



/(c) 



«7(i-" 2 ) 



2)3/2 



,-Pc 



(52) 



For higher dimensions (d > 2) we only quote the asymp- 
totic solution, 



/(c) 



c l~d/2 e -0c 



(53) 



which may also be obtained directly from (|47[) by ne- 
glecting the inhomogeneity, i.e. the gain term ^g am 

~ exp[-c 2 /a 2 ]. 

Comparison with the results of 32J for the WN-driven 
soft sphere models shows that the large-c behavior is ex- 
actly the same as that of the Maxwell model (y = 0), 
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FIG. 12: BGK with WN driving, in one and two dimensions, 
showing the exponentially decreasing tails. 



rich behavior of the Boltzmann equation, in particular it 
fails for hard interactions [y > 0). 





E=c 



FIG. 13: Isobestic points as either 6 is changed at constant v 
and a, or as v is changed at constant 8 and a. Top: d = 1; 
Bottom d = 2 (in 2 dimensions we plot the distribution of the 
energy E — c 2 ) . 



but the scaling solutions, displayed in Fig. fTJl are inde- 
pendent of v (since Eq. (|4"T)) is itself independent of ^), 
whether the scaling solution is a stable attracting state 
of a hard scattering model, or an unstable repelling state 
state of a soft scattering model. It shows therefore again 
that the BGK model is inadequate to model hard inter- 
actions. 

In summary, the simple linear BGK model, although 
displaying interesting features, such as power-law veloc- 
ity distribution tails, is far from being able to capture the 



VI. CONCLUDING REMARKS 

Within the framework of the nonlinear Boltzmann 
equation coupled to stochastic or deterministic driving 
forces and 'heat' baths, we have studied a general class of 
inelastic soft sphere models. Our approach encompasses 
a broad class of previously introduced models, from hard 
scatterers like inelastic hard spheres (and even very hard 
spheres [HI), to soft scatterers like Maxwell molecules, 
and even softer ones with v < 0, where v governs the 
dependence of collision frequency on relative velocity g 
through a term g v . 

We have shown that the velocity distribution /(c) has 
a stretched exponential tail cx exp(— c b ), when the non- 
equilibrium steady state is an attractive fixed point of 
the dynamics. In certain regions of model parameters 
{v, a, 9) where a denotes the restitution coefficient and 
9 is a friction parameter, we have reported important 
sub-leading corrections, where /(c) is found to be of the 
form c x exp(— f3v + (3'v b ). The comparison with high- 
precision numerical solutions of the Boltzmann equa- 
tion, obtained through Monte Carlo simulations (DSMC 
scheme), shows that neglecting these sub- leading correc- 
tions in a fitting procedure can lead to erroneous esti- 
mates of (3. Algebraic distributions emerge in cases of 
marginal stability (b — 0), and we have calculated the 
corresponding power law exponents. The high accuracy 
of our DSMC simulations have enabled us to verify the 
theoretical predictions for a wide range of parameter val- 
ues. 

The models studied here are partially amenable to an- 
alytical progress, but some features resist understanding. 
We conclude here by reporting one such feature, that is 
illustrated in Fig. [T3] We observe that all steady state 
rescaled velocity distributions, at fixed v and varying 9, 
pass through a common point. A similar property seems 
to hold when 9 is held fixed, and varying v. Such points, 
that can be coined "isobestic" , have already been ob- 
served in a different context (see e.g. section IV-E in 
reference [ll[), where their occurrence could not be ra- 
tionalized. 

Acknowledgements We would like to thank a referee 
for valuable suggestions. 
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